#Step1: Loadning relevant libraries
library(sf)
library(stars)
library(ggplot2)
sf_use_s2(FALSE)
#Step2: Read country borders
borders <- st_read(dsn="./data/boundaries/all-country-boundaries.shp", quiet = TRUE)
#Step3: Simplifying border to make them easier to use on our machine
borders<-st_simplify(borders, dTolerance=0.1)
#Step4: Creating a folder called "F101992"
target_dir <- "./data/luminosity/F101992"
#Step5: Create the target directory if it doesn't exist
if (!dir.exists(target_dir)) {
dir.create(target_dir)
}
#Step6: Unzip the file
untar("./data/luminosity/F101992.v4.tar", exdir = target_dir)
#Step7: Unzip the file further
R.utils::gunzip("./data/luminosity/F101992/F101992.v4b_web.avg_vis.tif.gz", remove = FALSE, overwrite=TRUE)
#Step8: Read the raster
f1992<-read_stars("./data/luminosity/F101992/F101992.v4b_web.avg_vis.tif")
#Step9: Make sure the two files have the same CRS
st_crs(f1992)<-st_crs(borders)
#Step10: Reduce the resolution of the raster: Higher dx/day values mean lower resolution
f1992b<-st_warp(f1992, st_as_stars(st_bbox(f1992), dx = 0.3, dy = 0.3))
#Step11: Rename the raster
names(f1992b)<-"lights_1992"
#Step12:Defining the countries defining the regions of interest
norway<-subset(borders, SOVEREIGNT == "Norway")
max_lat_y_eu<-st_bbox(norway)["ymax"]
greece<-subset(borders, SOVEREIGNT == "Greece")
min_lat_y_eu<-st_bbox(greece)["ymin"]
ukraine<-subset(borders, SOVEREIGNT == "Ukraine")
max_lon_x_eu<-st_bbox(ukraine)["xmax"]
portugal<-subset(borders, SOVEREIGNT == "Portugal")
min_lon_x_eu<-st_bbox(portugal)["xmin"]
#Step13: Making the map
fig1992<-ggplot()+
geom_stars(data=f1992b)+
scale_fill_gradientn(name = "", colours=c("black","white"))+
geom_sf(data=borders, linewidth = 0.1, fill = NA, color = "white", alpha=0.5)+
coord_sf(xlim = c(min_lon_x_eu-3, max_lon_x_eu+3), ylim = c(min_lat_y_eu-1, max_lat_y_eu-6))+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())L16: Rasters and Raster Operations 2
Zonal Statistics, Vectorizing, Cropping
Working with Luminosity 1992
Let us make a map of luminosity by focusing on Europe in 1992 and 2012
Working with Luminosity 2012
Working with Luminosity 2012
We can obviously put them side by side for comparison:
Germany
Let us first load the Germany county shape files.
Germany
Let us also load the Germany country shape file.
Germany
Let us now compare the two
Germany
We can now map the stars object
Germany
What if we tried to superimpose counties (kreise) on top of the raster?
Zonal Statiscs
This is for example what mean zonal statistics means in the context of luminosity
Zonal Statiscs
This is for example what mean zonal statistics maximum in the context of luminosity
Germany
We now look at the output:
Stars to Polygon
Let us now look at the output: the difference between starts and sf objects
Stars to Polygon
Visually, there is no difference between stars and sf
Germany: East vs. West
So, is there a difference between East and West?
Germany: East vs. West
So, is there a difference between East and West?
Step2: Subsetting based on that list
We can then do a subset just based on those states of NAME_1
Let us also plot them in the context of Germany.
Step3: Dissolving the county polygons to one polygon
We then dissolve the county polygons
This is what that union looks like:
Step4: Creating a 500m buffer to eliminate internal geometry errors
#Step5: Changing the coordinate system for Germany CRS
east_ger_diss2 = st_transform(east_ger_diss, 6875)
#Step6: Creating a 500m buffer to eliminate internal geometry errors
east_ger_diss3 <- st_buffer(east_ger_diss2, dist = 500)
#Step7: Turning back to WGS84
east_ger_diss4<-st_transform(st_as_sf(east_ger_diss3), st_crs(east_ger))Here is how the two polygons compare:
Step5: Create county centroids to identify their location in East or in the West
We then create the polygon centroid:
Step6: Perfoming the spatial join
We then add a column to indicate East and West.
Here is what the points look like after creating that column.
Step7: Merging the spatial attribute back to polygons
We then merge the region column from the points to the polygon.
#Step12: Dropping Geometry
ger_centr_joined_df<-st_drop_geometry(ger_centr_joined)
#Step13: Keeping only relevant variable resulting from the spatial merge
ger_centr_joined_df2<-subset(ger_centr_joined_df, select=c("NAME_2", "region"))
#Step14: Merging it back to the polygon
ger2c<-left_join(agg_sf2, ger_centr_joined_df2, by = c("NAME_2" = "NAME_2"))Here is what the output looks like:
Step 2: Turning the polygong to line
Let us turn the polygon into line
This is what the output looks like:
Step 3: Making neg. buffer to make the polygon size smaller so that we crop the internal eastern line
Here is the difference:
Step 4: Intersecting line with the smaller Germany polygon
We then intersect the line with the polygon:
Step 4: Intersecting line with the smaller Germany polygon
It seems that we have two lines:
Here is what they look like:
Step 4: Intersecting line with the smaller Germany polygon
Plotting them together
We will thus keep only line1.
Step 5: Calculating distance to line
Here is what the distances look like:
Comparison 1992
Let us draw a boxplot to see the differences:
Show the code
library(classInt)
# Desaturation libaries
library(colorspace)
library(colorblindr)
library(gridExtra)
library(viridis)
clint <- classIntervals(ger2d$lights_1992, style = "box")$brks
desaturation <- 1
# Get the original Viridis palette
original_palette <- viridis_pal()(100)
desaturated_palette <- desaturate(original_palette, desaturation)
# Plot with desaturated Viridis color scale
picx1_des<-ggplot() +
geom_sf(data = ger2d, aes(fill = lights_1992), color = "white", linewidth = 0.01) +
scale_fill_gradientn(colors = desaturated_palette, breaks = clint) +
guides(fill = guide_coloursteps(even.steps = FALSE,
show.limits = FALSE,
title = NULL,
barheight = unit(0.9, "npc"), # Relative to plot height
barwidth = unit(0.05, "npc"), # Relative to plot width
label.position = "left"))+
theme_bw()+
geom_sf(data = east_ger_diss4, color="red", linewidth=1, fill=NA)
picx2<-ggplot(ger2d, aes(x = region,
y = lights_1992,
color=region)) +
geom_boxplot()+
scale_y_continuous(breaks = clint,
labels = round(clint,2)) +
ylab(NULL) +
xlab(NULL) +
theme(axis.text = element_text(colour = "grey"),
panel.background = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank())
grid.arrange(picx1_des, picx2, ncol=2)
Comparison 2012
What about 2012?
Show the code
library(classInt)
# Desaturation libaries
library(colorspace)
library(colorblindr)
library(gridExtra)
clint <- classIntervals(ger2d$lights_2012, style = "box")$brks
desaturation <- 1
# Get the original Viridis palette
original_palette <- viridis_pal()(100)
desaturated_palette <- desaturate(original_palette, desaturation)
# Plot with desaturated Viridis color scale
picx1_des<-ggplot() +
geom_sf(data = ger2d, aes(fill = lights_2012), color = "white", linewidth = 0.01) +
scale_fill_gradientn(colors = desaturated_palette, breaks = round(clint,2)) +
guides(fill = guide_coloursteps(even.steps = FALSE,
show.limits = FALSE,
title = NULL,
barheight = unit(0.9, "npc"), # Relative to plot height
barwidth = unit(0.05, "npc"), # Relative to plot width
label.position = "left"))+
theme_bw()+
geom_sf(data = east_ger_diss4, color="red", linewidth=1, fill=NA)
picx2<-ggplot(ger2d, aes(x = region,
y = lights_2012,
color=region)) +
geom_boxplot()+
scale_y_continuous(breaks = clint,
labels = round(clint,2)) +
ylab(NULL) +
xlab(NULL) +
theme(axis.text = element_text(colour = "grey"),
panel.background = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank())
grid.arrange(picx1_des, picx2, ncol=2)
Difference
Let’s calculate the difference between 2012 and 1992.
Show the code
Comparison Difference
We can finally examine the change between 2012 and 1992.
Show the code
library(classInt)
# Desaturation libaries
library(colorspace)
library(colorblindr)
library(gridExtra)
clint <- classIntervals(ger2d$diff, style = "box")$brks
desaturation <- 0.2
# Get the original Viridis palette
original_palette <- viridis_pal()(100)
desaturated_palette <- desaturate(original_palette, desaturation)
# Plot with desaturated Viridis color scale
picx1_des<-ggplot() +
geom_sf(data = ger2d, aes(fill = diff), color = "white", linewidth = 0.01) +
scale_fill_viridis_c(option = "turbo", # You can choose other options like "A", "B", or "D"
#direction = -1, # Use -1 for reverse color order
breaks = clint,
guide = guide_coloursteps(even.steps = FALSE,
show.limits = FALSE,
title = NULL,
barheight = unit(0.9, "npc"), # Relative to plot height
barwidth = unit(0.05, "npc"), # Relative to plot width
label.position = "left"))+
theme_bw()+
geom_sf(data = line1, color="red", linewidth=2)
picx2<-ggplot(ger2d, aes(x = region,
y = diff,
color=region)) +
geom_boxplot()+
scale_y_continuous(breaks = clint,
labels = round(clint,2)) +
ylab(NULL) +
xlab(NULL) +
theme(axis.text = element_text(colour = "grey"),
panel.background = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank())
grid.arrange(picx1_des, picx2, ncol=2)
Refining the Analysis
One criticism to our analysis is that we are comparing units that are too different.
For example, it feels strange to compare units like Lörrach or Uckermark.
Show the code
ger_sel<-subset(ger2d, NAME_2=="Lörrach" | NAME_2=="Uckermark")
ggplot() +
geom_sf(data = ger0) +
geom_sf(data = ger2d, aes(fill = diff), size = 0) +
geom_sf(data = east_ger_diss4, fill = NA, color = "red", linewidth = 1) +
scale_fill_viridis(option="turbo",
name = "Lum. 2012-\nLum. 1992",
guide = guide_coloursteps(even.steps = FALSE))+
geom_sf_label(data=ger_sel, aes(label=NAME_2), size=1.5)
Refining the Analysis
Counties like Lörrach or Uckermark may be different for other reasons that have nothing to do with exposure to communism.
Thus, it may be more appropriate to examine counties that are close the East-West border
Show the code
ggplot()+
geom_sf(data=ger0, fill=NA)+
geom_sf(data=ger_100km, aes(fill=lights_2012-lights_1992))+
scale_fill_viridis(option="turbo",
name = "Lum. 2012-\nLum. 1992",
guide = guide_coloursteps(even.steps = FALSE))+
geom_sf(data=line1, color="red", linewidth=1.5)+
coord_sf(ylim=c(ymin,ymax),
xlim=c(xmin, xmax))
Comparison 1992
We can now perform the analysis only based on the 100km sample around the border
Show the code
library(classInt)
# Desaturation libaries
library(colorspace)
library(colorblindr)
library(gridExtra)
clint <- classIntervals(ger_100km$lights_1992, style = "box")$brks
desaturation <- 1
# Get the original Viridis palette
original_palette <- viridis_pal()(100)
desaturated_palette <- desaturate(original_palette, desaturation)
# Plot with desaturated Viridis color scale
picx1_des<-ggplot() +
geom_sf(data = ger_100km, aes(fill = lights_1992), color = "white", linewidth = 0.01) +
scale_fill_gradientn(colors = desaturated_palette, breaks = clint) +
guides(fill = guide_coloursteps(even.steps = FALSE,
show.limits = FALSE,
title = NULL,
barheight = unit(0.9, "npc"), # Relative to plot height
barwidth = unit(0.05, "npc"), # Relative to plot width
label.position = "left"))+
theme_bw()+
geom_sf(data = line1, color="red", linewidth=2)
picx2<-ggplot(ger_100km, aes(x = region,
y = lights_1992,
color=region)) +
geom_boxplot()+
scale_y_continuous(breaks = clint,
labels = round(clint,2)) +
ylab(NULL) +
xlab(NULL) +
theme(axis.text = element_text(colour = "grey"),
panel.background = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank())
grid.arrange(picx1_des, picx2, ncol=2)
Comparison 2012
We can also perform the analysis only based on the 100km sample around the border for 2012
Show the code
library(classInt)
# Desaturation libaries
library(colorspace)
library(colorblindr)
library(gridExtra)
clint <- classIntervals(ger_100km$lights_2012, style = "box")$brks
desaturation <- 1
# Get the original Viridis palette
original_palette <- viridis_pal()(100)
desaturated_palette <- desaturate(original_palette, desaturation)
# Plot with desaturated Viridis color scale
picx1_des<-ggplot() +
geom_sf(data = ger_100km, aes(fill = lights_2012), color = "white", linewidth = 0.01) +
scale_fill_gradientn(colors = desaturated_palette, breaks = clint) +
guides(fill = guide_coloursteps(even.steps = FALSE,
show.limits = FALSE,
title = NULL,
barheight = unit(0.9, "npc"), # Relative to plot height
barwidth = unit(0.05, "npc"), # Relative to plot width
label.position = "left"))+
theme_bw()+
geom_sf(data = line1, color="red", linewidth=2)
picx2<-ggplot(ger_100km, aes(x = region,
y = lights_2012,
color=region)) +
geom_boxplot()+
scale_y_continuous(breaks = clint,
labels = round(clint,2)) +
ylab(NULL) +
xlab(NULL) +
theme(axis.text = element_text(colour = "grey"),
panel.background = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank())
grid.arrange(picx1_des, picx2, ncol=2)
Comparison Difference
Show the code
library(classInt)
# Desaturation libaries
library(colorspace)
library(colorblindr)
library(gridExtra)
ger_100km$diff<-ger_100km$lights_2012-ger_100km$lights_1992
clint <- classIntervals(ger_100km$diff, style = "box")$brks
desaturation <- 0.2
# Get the original Viridis palette
original_palette <- viridis_pal()(100)
desaturated_palette <- desaturate(original_palette, desaturation)
# Plot with desaturated Viridis color scale
picx1_des<-ggplot() +
geom_sf(data = ger_100km, aes(fill = diff), color = "white", linewidth = 0.01) +
scale_fill_viridis_c(option = "turbo", # You can choose other options like "A", "B", or "D"
#direction = -1, # Use -1 for reverse color order
breaks = clint,
guide = guide_coloursteps(even.steps = FALSE,
show.limits = FALSE,
title = NULL,
barheight = unit(0.9, "npc"), # Relative to plot height
barwidth = unit(0.05, "npc"), # Relative to plot width
label.position = "left"))+
theme_bw()+
geom_sf(data = line1, color="red", linewidth=2)
picx2<-ggplot(ger_100km, aes(x = region,
y = diff,
color=region)) +
geom_boxplot()+
scale_y_continuous(breaks = clint,
labels = round(clint,2)) +
ylab(NULL) +
xlab(NULL) +
theme(axis.text = element_text(colour = "grey"),
panel.background = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank())
grid.arrange(picx1_des, picx2, ncol=2)